Introduction

This is the code to reproduce Figure 3S Additional longitudinal and repertoire analysis of T cells across timepoints.: Fig3S C-H. To obtain the data object used in this notebook, please run 01_TCR_Data_Analysis.Rmd.

Package loading

library(gdata)
library(Seurat)
library(reshape2)
library(ggplot2)
library(ggpubr)
library(cowplot)
library(scRepertoire)
library(circlize)
library(Startrac)
library(stringr)
library(plyr)
library(dplyr)
library(ggh4x)
library(ggbeeswarm)
library(monocle)

Loading colors

cols_patient =c("Patient 1"= "aquamarine", "Patient 2"="lightpink", "Patient 3"="yellow1", "Patient 4"="skyblue", "Patient 5"="sienna1")
cols_patient2 <- c("aquamarine", "#6FDFBA", "lightpink", "#DF9FA9", "yellow1", "#DFDC00", "skyblue", "#76B4CF","sienna1", "#DF733E")
cols_timepoint<- c("IP"="#4E6AAB","Peak"="#e78ac3")
cols_anno <- c("CD4+ Naive T cells"= "#33A02C",
                     "CD4+ CEntral/Effector memory T cells (CM/EM)"="#B2DF8A",
                     "CD8+ cytotoxic T cells"="#185B88",
                     "CD8+ Effector T cells (E)"="#1F78B4",
                     "CD8+ Eff/Mem T cells (EM)"="#A6CEE3",
                     "Early prolif: MCM3/5/7+ PCNA+ T cells"="#FB9A99",
                     "Late prolif: histones enriched MKI67+ T cells"="#E31A1C",
                     "Late prolif: CCNB1/2+ CDK1+ T cells"="#CAB2D6", 
                     "Late prolif: STMN1+ BIRC5+"="#FDBF6F",
                     "Ribosomal/Mitochondrial/Degraded cells"="#FF7F00",
                     "gamma-delta T cells"="#6A3D9A")
strip <- strip_themed(background_x = elem_list_rect(fill = c("#A6CEE3", "#1F78B4","#185B88", "#B2DF8A","#33A02C", "#CAB2D6", "#FB9A99", "#E31A1C", "#FDBF6F", "#FF7F00" )))
col_clono <- c("Hyperexpanded (100 < X <= 500)"="#810F7C", "Large (20 < X <= 100)"="#8856A7","Medium (5 < X <= 20)"= "#8C96C6","Small (1 < X <= 5)"= "#9EBCDA", "Single (0 < X <= 1)"="#BFD3E6", "No clonotype detected"="gray69")

Data import

combined_infusion <- readRDS("../data/combined_files/combined_infusion_CARPOS.rds")
combined_infusion<- addVariable(combined_infusion, variable.name  = "patient", 
                variables = c("patient1", "patient2","patient3", "patient4","patient5"))
table_combined_infusion<- clonalDiversity(combined_infusion, 
                cloneCall = "gene", 
                group.by = "sample", 
                n.boots = 100, exportTable = T) #diversity_infusion
combined_peak <- readRDS("../data/combined_files/combined_peak_CARPOS.rds")
combined_peak<- addVariable(combined_peak, variable.name = "patient", 
                variables = c("patient1", "patient2","patient3", "patient4","patient5"))
table_combined_peak<- clonalDiversity(combined_peak, 
                cloneCall = "gene", 
                group.by = "sample", 
                n.boots = 100, exportTable = T) #diversity_peak
table_combined_infusion$timepoint <- "Infusion"
table_combined_peak$timepoint <- "Peak"

table_combined_infusion$patient <- c("Patient 1", "Patient 2", "Patient 3", "Patient 4", "Patient 5")
table_combined_peak$patient <- c("Patient 1", "Patient 2", "Patient 3", "Patient 4", "Patient 5")

table_combined <- rbind(table_combined_infusion,table_combined_peak)
table_combined$sample <- factor(table_combined$sample)
seurat <- readRDS("../data/MENENDEZ_DEF.rds")
seurat_carneg <- subset(seurat, subset = Class1 == "CAR-")
seurat_carpos <- subset(seurat, subset = Class1 == "CAR+")

seurat_carneg <- subset(seurat_carneg, subset = cloneSize != "No clonotype detected")
seurat_carneg <- subset(seurat_carneg, subset = cloneSize != "Single (0 < X <= 1)")
seurat_carneg_nogd <- subset(seurat_carneg, subset = annotation != "gamma-delta T cells")

seurat_carpos <- subset(seurat_carpos, subset = cloneSize != "No clonotype detected")
seurat_carpos <- subset(seurat_carpos, subset = cloneSize != "Single (0 < X <= 1)")
seurat_carpos_nogd <- subset(seurat_carpos, subset = annotation != "gamma-delta T cells")

seurat_carneg_IP <- subset(seurat_carneg_nogd, subset = Timepoint == "IP")
seurat_carneg_PEAK <- subset(seurat_carneg_nogd, subset = Timepoint == "Peak")
data_pt1 <- subset(seurat_carpos_nogd, subset = Patient_id == "patient1")
data_pt1_I <- subset(data_pt1, subset = Timepoint == "IP")
data_pt1_P <- subset(data_pt1, subset = Timepoint == "Peak")

data_pt2 <- subset(seurat_carpos_nogd, subset = Patient_id == "patient2")
data_pt2_I <- subset(data_pt2, subset = Timepoint == "IP")
data_pt2_P <- subset(data_pt2, subset = Timepoint == "Peak")

data_pt3 <- subset(seurat_carpos_nogd, subset = Patient_id == "patient3")
data_pt3_I <- subset(data_pt3, subset = Timepoint == "IP")
data_pt3_P <- subset(data_pt3, subset = Timepoint == "Peak")

data_pt4 <- subset(seurat_carpos_nogd, subset = Patient_id == "patient4")
data_pt4_I <- subset(data_pt4, subset = Timepoint == "IP")
data_pt4_P <- subset(data_pt4, subset = Timepoint == "Peak")

data_pt5 <- subset(seurat_carpos_nogd, subset = Patient_id == "patient5")
data_pt5_I <- subset(data_pt5, subset = Timepoint == "IP")
data_pt5_P <- subset(data_pt5, subset = Timepoint == "Peak")
data2<- readRDS("../data/dataordered_V2.rds")

Plot Fig3S C

occRep_ip<-clonalOccupy(seurat_carneg_IP, x.axis = "annotation", exportTable =T)
occRep_peak <-clonalOccupy(seurat_carneg_PEAK, x.axis = "annotation", exportTable =T)

occRep_ip$Timepoint <- "IP"
occRep_peak$Timepoint <- "PEAK"

occRep <- rbind(occRep_ip, occRep_peak)
occRep$annotation <- factor(occRep$annotation, levels = 
                                    c("CD8+ Eff/Mem T cells (EM)", 
                                      "CD8+ Effector T cells (E)", 
                                      "CD8+ cytotoxic T cells", 
                                      "CD4+ CEntral/Effector memory T cells (CM/EM)", 
                                      "CD4+ Naive T cells", "gamma-delta T cells", 
                                      "Late prolif: CCNB1/2+ CDK1+ T cells", 
                                      "Early prolif: MCM3/5/7+ PCNA+ T cells", 
                                      "Late prolif: histones enriched MKI67+ T cells", 
                                      "Late prolif: STMN1+ BIRC5+",
                                      "Ribosomal/Mitochondrial/Degraded cells"))
occRep$cloneSize <- factor(occRep$cloneSize, levels = c("Small (1 < X <= 5)", "Medium (5 < X <= 20)", "Large (20 < X <= 100)", "Hyperexpanded (100 < X <= 500)"))


 # Draw barplot with grouping & stacking
ggplot(occRep, aes(x = Timepoint, y = n, fill = cloneSize)) + 
  geom_bar(stat = "identity", position = "stack") +
  facet_grid(~ annotation) +
  scale_fill_manual(values = col_clono) + facet_grid2(~ annotation, strip = strip) +
  theme(strip.text = element_text(colour = NA),
        axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=0.5))

Plot Fig3S D

I <- getCirclize(data_pt1_I, group.by = "annotation")
I<-chordDiagram(I, self.link = 1, grid.col =cols_anno , annotationTrack = c("grid", "axis"))
title("IP Patient 1")

P <- getCirclize(data_pt1_P, group.by = "annotation")
P<-chordDiagram(P, self.link = 1, grid.col =cols_anno, annotationTrack = c("grid", "axis") )
title("Peak Patient 1")

I <- getCirclize(data_pt2_I, group.by = "annotation")
I<-chordDiagram(I, self.link = 1, grid.col =cols_anno , annotationTrack = c("grid", "axis"))
title("IP Patient 2")

P <- getCirclize(data_pt2_P, group.by = "annotation")
P<-chordDiagram(P, self.link = 1, grid.col =cols_anno, annotationTrack = c("grid", "axis") )
title("Peak Patient 2")

I <- getCirclize(data_pt3_I, group.by = "annotation")
I<-chordDiagram(I, self.link = 1, grid.col =cols_anno , annotationTrack = c("grid", "axis"))
title("IP Patient 3")

P <- getCirclize(data_pt3_P, group.by = "annotation")
P<-chordDiagram(P, self.link = 1, grid.col =cols_anno, annotationTrack = c("grid", "axis") )
title("Peak Patient 3")

I <- getCirclize(data_pt4_I, group.by = "annotation")
I<-chordDiagram(I, self.link = 1, grid.col =cols_anno , annotationTrack = c("grid", "axis"))
title("IP Patient 4")

P <- getCirclize(data_pt4_P, group.by = "annotation")
P<-chordDiagram(P, self.link = 1, grid.col =cols_anno, annotationTrack = c("grid", "axis") )
title("Peak Patient 4")

I <- getCirclize(data_pt5_I, group.by = "annotation")
I<-chordDiagram(I, self.link = 1, grid.col =cols_anno , annotationTrack = c("grid", "axis"))
title("IP Patient 5")

P <- getCirclize(data_pt5_P, group.by = "annotation")
P<-chordDiagram(P, self.link = 1, grid.col =cols_anno, annotationTrack = c("grid", "axis") )
title("Peak Patient 5")

Plot Fig3S E

text <- element_text(color = "black", size = 16)
text.lab <- element_text(color = "black", size = 12)
text.lab2 <- element_text(color = "black", size = 10)
text.lab3 <- element_text(color = "black", size = 9)

aux_df <- data.frame(barcodes=colnames(data2),idents=data2$annotation,Pseudotime=data2$Pseudotime, Timepoint = data2$Timepoint)

ident_trajectory <- ggplot(aux_df,aes(x=idents,y=Pseudotime,col=Timepoint)) +
  geom_quasirandom(alpha=.8) + coord_flip() + theme_classic() +
  scale_color_manual(values=cols_timepoint) + theme(text = element_text(size = 14)) 

ident_trajectory

Plot Fig3S F

ordering.genes <- VariableFeatures(seurat)
marker_genes <- row.names(subset(fData(data2),
                   gene_short_name %in% ordering.genes))
data2 <- reduceDimension(data2,
                         max_components = 2,
                         norm_method = 'log',
                         num_dim = 3,
                         reduction_method = 'tSNE',
                         verbose = T)
data2 <- clusterCells(data2, verbose = F, num_clusters = 10)
diff_test_res <- differentialGeneTest(data2[marker_genes,],
                                      fullModelFormulaStr = "~sm.ns(Pseudotime)",
                                      cores=4)
sig_gene_names <- row.names(subset(diff_test_res, qval < 0.000000000000000000000000000000000000001))
sig_gene_names <- unique(c("STMN1", "CDC20", "HIST2H2BF", "UBE2C", "GNLY", "GZMK", "CCNB1", "MKI67", "KLRG1", "KLRB1", "CD69", "GZMH", "HMGB1", "PCNA","GZMB", "IFITM1", "CD27", "LAG3", "MKI67", "TOP2A", "CDK1", "CXCR3", "GZMK", "NKG7", "HSPA1A", "HSPA1B", "PRF1", "TRDV1", "TRGV2" ))
plot_pseudotime_heatmap(data2[sig_gene_names,],
                num_clusters = 3,
                cores = 1,
                show_rownames = T,
                hmcols = colorRampPalette(colors = c("blue", "white", "red"))(62))

Plot Fig3S G

cols2=c("CD4+ CEntral/Effector memory T cells (CM/EM)"="white",
               "CD4+ Naive T cells"="white",
               "CD8+ cytotoxic T cells"="#185B88",
               "CD8+ Eff/Mem T cells (EM)"="white",
               "CD8+ Effector T cells (E)"="white",
               "Early prolif: MCM3/5/7+ PCNA+ T cells"="white",
               "Late prolif: CCNB1/2+ CDK1+ T cells"="white",
               "Late prolif: histones enriched MKI67+ T cells"="white",
               "Late prolif: STMN1+ BIRC5+"="white",
               "Ribosomal/Mitochondrial/Degraded cells"="white")
data_pt1 <- subset(seurat_carneg_nogd, subset = Patient_id == "patient1")
data_pt2 <- subset(seurat_carneg_nogd, subset = Patient_id == "patient2")
data_pt3 <- subset(seurat_carneg_nogd, subset = Patient_id == "patient3")
data_pt4 <- subset(seurat_carneg_nogd, subset = Patient_id == "patient4")
data_pt5 <- subset(seurat_carneg_nogd, subset = Patient_id == "patient5")


pt1 <-alluvialClones(data_pt1, cloneCall = "CTstrict", 
                   y.axes = c("annotation", "cloneType", "Timepoint"), 
                   color = "annotation") + scale_fill_manual(values=cols2)+theme(legend.position = "none") +
  scale_y_continuous(limits = c(0, 1650))
pt2 <-alluvialClones(data_pt2, cloneCall = "CTstrict", 
                   y.axes = c("annotation", "cloneType", "Timepoint"), 
                   color = "annotation") + scale_fill_manual(values=cols2)+theme(legend.position = "none") +
  scale_y_continuous(limits = c(0, 1650))
pt3 <-alluvialClones(data_pt3, cloneCall = "CTstrict", 
                   y.axes = c("annotation", "cloneType", "Timepoint"), 
                   color = "annotation") + scale_fill_manual(values=cols2)+theme(legend.position = "none") +
  scale_y_continuous(limits = c(0, 1650))
pt4 <-alluvialClones(data_pt4, cloneCall = "CTstrict", 
                   y.axes = c("annotation", "cloneType", "Timepoint"), 
                   color = "annotation") + scale_fill_manual(values=cols2)+theme(legend.position = "none") +
  scale_y_continuous(limits = c(0, 1650))
pt5 <-alluvialClones(data_pt5, cloneCall = "CTstrict", 
                   y.axes = c("annotation", "cloneType", "Timepoint"), 
                   color = "annotation") + scale_fill_manual(values=cols2)+theme(legend.position = "none") + scale_y_continuous(limits = c(0, 1650))

ggarrange(pt1 + facet_grid(. ~"Patient 1"), 
          pt2 + 
               theme(axis.text.y = element_blank(),
                     axis.ticks.y = element_blank(),
                     axis.title.y = element_blank(),
                     plot.margin = margin(r = 1, l = 1) )+ facet_grid(. ~"Patient 2"),
          pt3 + 
               theme(axis.text.y = element_blank(),
                     axis.ticks.y = element_blank(),
                     axis.title.y = element_blank(),
                     plot.margin = margin(r = 1, l = 1) )+ facet_grid(. ~"Patient 3"), 
          pt4 + 
               theme(axis.text.y = element_blank(),
                     axis.ticks.y = element_blank(),
                     axis.title.y = element_blank(),
                     plot.margin = margin(r = 1, l = 1) )+ facet_grid(. ~"Patient 4"), 
          pt5 + 
               theme(axis.text.y = element_blank(),
                     axis.ticks.y = element_blank(),
                     axis.title.y = element_blank(),
                     plot.margin = margin(r = 1, l = 1) )+ facet_grid(. ~"Patient 5"),  align = "hv", ncol=5)

Plot Fig3S H

seurat_carpos <- subset(seurat, subset = Class1 == "CAR+")

table_unique_clono <- scRepertoire:::.expression2List(seurat_carpos, split.by = "Sample_id")
table_unique_clono <-clonalQuant(table_unique_clono,  exportTable = T)
table_unique_clono$patient <- revalue(table_unique_clono$values, c("patient1_IP"="Patient 1", "patient1_Peak"="Patient 1", "patient2_IP"="Patient 2", "patient2_Peak"="Patient 2", "patient3_IP"="Patient 3", "patient3_Peak"="Patient 3", "patient4_IP"="Patient 4", "patient4_Peak"="Patient 4", "patient5_IP"="Patient 5", "patient5_Peak"="Patient 5" ))
table_unique_clono$patient <- factor(table_unique_clono$patient)
table_unique_clono$timepoint <- rep(c("IP","Peak"),5)

ggplot(table_unique_clono, aes(fill=values, y=contigs, x=patient)) + 
    geom_bar(position="dodge", stat="identity")+ scale_fill_manual(values=cols_patient2) +
  geom_text(aes(label=timepoint), position = position_dodge(width = .9), vjust = -0.3, color = "black") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))

Session Info

sessionInfo()
## R version 4.2.0 (2022-04-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur/Monterey 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] splines   stats4    stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] monocle_2.24.0      DDRTree_0.1.5       irlba_2.3.5.1       VGAM_1.1-9          Biobase_2.58.0      BiocGenerics_0.44.0 Matrix_1.6-5        ggbeeswarm_0.7.2    ggh4x_0.2.8         dplyr_1.1.4         plyr_1.8.9          stringr_1.5.1       Startrac_0.1.0      circlize_0.4.16     scRepertoire_2.0.4  cowplot_1.1.3       ggpubr_0.6.0        ggplot2_3.5.1       reshape2_1.4.4      Seurat_5.0.1        SeuratObject_5.0.1  sp_2.1-2            gdata_3.0.0        
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.4                  spatstat.explore_3.2-5      reticulate_1.34.0           tidyselect_1.2.1            htmlwidgets_1.6.4           docopt_0.7.1                combinat_0.0-8              grid_4.2.0                  Rtsne_0.17                  munsell_0.5.1               codetools_0.2-20            ica_1.0-3                   future_1.33.2               miniUI_0.1.1.1              withr_3.0.0                 fastICA_1.2-4               spatstat.random_3.2-2       colorspace_2.1-0            progressr_0.14.0            highr_0.10                  knitr_1.41                  ggalluvial_0.12.5           rstudioapi_0.16.0           SingleCellExperiment_1.20.1 ROCR_1.0-11                 ggsignif_0.6.4              tensor_1.5                  listenv_0.9.1               labeling_0.4.3              MatrixGenerics_1.10.0       slam_0.1-50                 GenomeInfoDbData_1.2.9      polyclip_1.10-6             pheatmap_1.0.12             farver_2.1.1                parallelly_1.36.0           vctrs_0.6.5                 generics_0.1.3              xfun_0.41                   doParallel_1.0.17           R6_2.5.1                    GenomeInfoDb_1.34.9        
##  [43] clue_0.3-65                 graphlayouts_1.0.1          bitops_1.0-7                spatstat.utils_3.0-4        cachem_1.0.8                DelayedArray_0.24.0         promises_1.2.1              scales_1.3.0                ggraph_2.1.0                beeswarm_0.4.0              gtable_0.3.5                globals_0.16.3              goftest_1.2-3               spam_2.10-0                 tidygraph_1.3.0             MatrixModels_0.5-3          rlang_1.1.2                 GlobalOptions_0.1.2         rstatix_0.7.2               lazyeval_0.2.2              spatstat.geom_3.2-7         broom_1.0.5                 yaml_2.3.8                  abind_1.4-5                 backports_1.4.1             httpuv_1.6.13               tools_4.2.0                 cubature_2.1.0              jquerylib_0.1.4             RColorBrewer_1.1-3          ggdendro_0.2.0              ggridges_0.5.6              hash_2.2.6.3                Rcpp_1.0.11                 zlibbioc_1.44.0             purrr_1.0.2                 RCurl_1.98-1.13             deldir_2.0-2                GetoptLong_1.0.5            pbapply_1.7-2               viridis_0.6.5               S4Vectors_0.36.2           
##  [85] zoo_1.8-12                  iNEXT_3.0.1                 SummarizedExperiment_1.28.0 ggrepel_0.9.4               cluster_2.1.6               magrittr_2.0.3              data.table_1.15.4           RSpectra_0.16-1             scattermore_1.2             SparseM_1.81                lmtest_0.9-40               RANN_2.6.1                  truncdist_1.0-2             fitdistrplus_1.1-11         matrixStats_1.2.0           gsl_2.1-8                   patchwork_1.2.0.9000        mime_0.12                   evaluate_0.23               xtable_1.8-4                sparsesvd_0.2-2             shape_1.4.6.1               fastDummies_1.7.3           IRanges_2.32.0              gridExtra_2.3               HSMMSingleCell_1.18.0       compiler_4.2.0              tibble_3.2.1                crayon_1.5.2                KernSmooth_2.23-22          htmltools_0.5.7             later_1.3.2                 tidyr_1.3.0                 tweenr_2.0.2                ComplexHeatmap_2.14.0       MASS_7.3-57                 leidenbase_0.1.27           car_3.1-2                   cli_3.6.2                   evd_2.3-6.1                 parallel_4.2.0              dotCall64_1.1-1            
## [127] igraph_1.5.1                GenomicRanges_1.50.2        pkgconfig_2.0.3             plotly_4.10.4               spatstat.sparse_3.0-3       foreach_1.5.2               vipor_0.4.7                 bslib_0.4.1                 stringdist_0.9.12           XVector_0.38.0              digest_0.6.33               sctransform_0.4.1           RcppAnnoy_0.0.21            spatstat.data_3.0-4         rmarkdown_2.18              leiden_0.4.3.1              uwot_0.1.16                 evmix_2.12                  shiny_1.8.1.1               gtools_3.9.5                quantreg_5.97               rjson_0.2.21                lifecycle_1.0.4             nlme_3.1-164                jsonlite_1.8.8              carData_3.0-5               limma_3.54.2                viridisLite_0.4.2           fansi_1.0.6                 pillar_1.9.0                lattice_0.22-5              fastmap_1.1.1               httr_1.4.7                  survival_3.5-7              glue_1.6.2                  qlcMatrix_0.9.7             iterators_1.0.14            png_0.1-8                   ggforce_0.4.1               stringi_1.8.3               sass_0.4.8                  RcppHNSW_0.5.0             
## [169] future.apply_1.11.2